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— The method of geodesic deviations provides analytic approximations to geodesies in arbitrary back- 

ground space-times. As such the method is a useful tool in many practical situations. In this note we 
point out some subtleties in the application of the method related to secular motions, in first as well 
as in higher order. In particular we work out the general second-order contribution to bound orbits in 
Schwarzschild space-time and show that it provides very good analytical results all the way up to the 
t-H innermost stable circular orbit. 
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1 Geodesic deviations 



According to General Relativity, in a fixed background space-time compact objects (test 
masses) move on geodesies. However, except for the simplest cases, their orbits can be cal- 
culated from the geodesic equation only in certain approximations. An often- used method is 
provided by the post-newtonian approximation scheme, which starts from the non-relativistic 
orbit, and then systematically calculates special and general relativistic corrections [USE]. 

A different approximation scheme is provided by the method of geodesic deviations [1, 2j. 
This is a manifestly covariant method, which can be extended to include other interactions, 
e.g. in Einstein-Maxwell theory [HE] and non-abelian backgrounds [6J, or the effects of spin 
[8]. Moreover, the method can be extended to arbitrary precision by taking into account 
higher-order deviations [91 [10] . In the literature cited [I] -[10], the method has been applied 
e.g. to calculate particle orbits in pp- waves, and Schwarzschild, Reissner-Nordstr0m and Kerr 
space-times. 

In this paper we consider the application of the geodesic deviation method to pure gravity, 
pointing out some subtleties which arise already at linear order, but which become of special 
relevance when extending the method to include higher-order corrections. Similar subtleties 
also arise in newtonian gravity, and have been adressed in newtonian perturbation theory 
already in the 19th century [TU [T2"] . 

The geodesic deviation method starts from a known reference geodesic, and then computes 
neighboring geodesies by determining the space-time vector connecting the points of the 
reference orbit with points on the unknown geodesic. Let x^^u] represent a continuous 
family of geodesies, the proper time r acting as affine parameter, and a labeling the geodesies 
in the family. Let x^(t) = x^[r; 0] be a known geodesic: 

D 2 x^ _d?x> 1 - ^dx* dx v _ 

where f = r[x(r)] represents the connection evaluated on the geodesic x[t]. Neighboring 
geodesies are then found from the expansion 

X ^[ T - a] = x»(t) + an^{r) + ~ ctW(t) + ... =^+(j/ + ^ a 2 (hP - r^nV) + ... (2) 

In this expansion the vectors and k^ are defined covariantly as 



da 



=o 



Da 



da 



+ r A >V. (3) 



The scale can be set, for example, by defining a to represent the proper distance along the 
curve x M [0; a]: 

da 2 = g fJ/l/ dx tl dx u \ T=0 => 9^n tl n u \ T=0 = 1, (4) 

where, as for the connection, = gf„^[x(r)] is the metric evaluated on the reference geodesic 
x^{t). Of course, as the distance between geodesies varies the normalization of n is in general 
not preserved as a function of proper time. Specifically, the change of the geodesic deviation 
vectors (n, k, ...) is determined by the geodesic deviation equations 

D 2 n n 



Dt 2 Kv^n" = 0, 



D 2 k fl / \ Dn K 



(5) 
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and higher-order generalizations [9], where is the four- velocity along the reference geodesic: 

(6) 



dx^- dx^ 
fir dr 



a=0 



By construction, these equations are successively linear in the perturbation vectors (n, k, ...); 
however, whereas the equation for the first-order perturbation is a homogeneous linear dif- 
ferential equation, the equation for the higher-order corrections are inhomogeneous linear 
differential equations. The inhomogeneous terms are polynomial in the lower-order perturba- 
tions, such that e.g. the second-order perturbation k is determined by terms quadratic in the 
first-order perturbation n. In these inhomogeneous quadratic terms n is to be considered as 
a known vector, having been solved from the homogeneous first-order equation. 

Another obvious property of these equations is, that the linear terms on the left-hand side 
are all of the same form [9]. This implies that the general solution of the equation for the 
second- and higher-order perturbations is some particular solution plus an arbitrary solution 
of the homogeneous equation, i.e. a first-order solution. 

Considering the first-order equation, we observe that its general form is that of a para- 
metric oscillator, with a proper-time dependent driving force linear in the amplitude n, the 
full proper-time dependence being determined by the curvature tensor and the four-velocity 
on the reference geodesic. The full power of the geodesic deviation method is developed if 
the manifold considered admits a reference geodesic for which the parametric force equation 
is solvable. Such solutions can be found in many important cases, such as the Schwarzschild, 
Reissner-Nordstr0m and Kerr solutions in four dimensional space-time. Below we consider 
in particular the Schwarzschild geometry, but the methodological observations hold quite 
generally. 



2 First-order perturbations 

The linear first-order equation is the key to all higher-order ones, as it provides both inho- 
mogenous terms for, and the general homogeneous solutions of, the higher-order equations. In 
spaces with symmetries there are some obvious solutions. Indeed, if £^(x) is a Killing vector 
field: 

Sw + Zm» = °> ( 7 ) 



it generates an isometry of the metric and it follows quite simply that it also satisfies the 
geodesic deviation equation: 

Rx„y x u K e, (8) 



D 2 ^ _ 

D T 2 ' -"A'".' 



along any geodesic. For example, in exterior Schwarzschild space-time, which is static and 
spherically symmetric, any geodesic is turned into another geodesic by a time translation 
or a rotation. By using this property one can in fact reduce the class of geodesies to be 
investigated by simply modding out part of the rotations, and considering only geodesies in 
the equatorial plane, as in the standard textbook treatment [DI2JEE3]. When discussing the 
geodesies of Schwarzschild space-time below we do the same. 

Another interesting situation occurs, if geodesies are simultaneously lines of Killing flow, 
i.e. if a Killing vector is transported parallel to itself along a geodesic. Then the tangent vector 
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(four- velocity) along such a geodesic is itself a Killing vector and its norm is constant along 
the geodesic 

^ 2 = 2a^ = -2a^ = 0' (9) 

An example is provided by circular equatorial orbits in Schwarzschild geometry (our conven- 
tions are summarized in appendix |A|) , which are characterized by a four- velocity 

« M = -eotf + 4$. (10) 

Here (£t, £^) are the Killing vectors generating time translations and axial rotations, and (e, £) 
are the corresponding constants of motion defined in appendix [Aj eq. (91), which on circular 
orbits take the values 



Explicitly, the Killing vectors for time translations and axial rotations are represented by the 
differential operators 

Zt = &d ll = g tt d t , ^ = ^ = g^8 v . (12) 



These Killing vectors are orthogonal: Ct ' C<p = 0) an d eqs. (11) imply that the linear combi- 
nation ( 10 ) is time- like and normalized: 

u 2 =el$+lie H , = -l- (13) 

Returning to the general discussion, if the tangent vector generates an isometry, the metric 
and connection are constant along the geodesic. With such a geodesic as a reference, the 
geodesic deviation equations then reduce to ordinary linear second-order differential equations 
with constant co-efficients. In particular, for bound orbits these equations have standard 
solutions with real eigenfrequencies [H [9] . 

To find explicit solutions, it is convenient to write out the first-order equation in its non- 
manifestly covariant form 

^ + 2n x F x ^ + u K u x d v t K £ n" = 0. (14) 

As for the special reference orbits the coefficients are constant, the generic solution for bound 
orbits will be periodic: 

n per = n c cos UJT + n s s i n w ' r i (15) 

with constant amplitudes (n c ,n s ). The eigenfrequencies uj are found by diagonalizing the 
characteristic matrix for the differential equations (14), defined by the set of linear equations 

(16) 



-u?n^ + 2cou x T x y s + u K u x d v F K £ < = 0, 



-uj 2 n% - 2uu x T x ^n u c + u K u x d v Y R £ < = 0. 

It follows, that in n space-time dimensions the characteristic equation in general has 2n roots. 
However, even if the reference orbit is strictly bound (i.e., it is enclosed in a finite region of 
space) and the eigenfrequencies are guaranteed to be real, there can still be zero-modes. These 
zero modes take the form of polynomials in t, rather than periodic functions of the type ( 15 ). 
Being non-periodic, they describe secular motions. The importance of, and how to deal with, 
such zero modes is addressed in the following. 
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3 First-order geodesies in Schwarzschild space-time 



Before continuing the general discussion, we analyze in more detail the example of orbits near 
a circular orbit in exterior Schwarzschild space-time. We have already seen that for such 
orbits we can reduce the problem from a four-dimensional one to a three-dimensional one by 
restriction to the equatorial plane. For motions in the equatorial plane, eq. (14) takes the 
form 



where 



/ 
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dr dr 
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(17) 



(18) 



Then the linear equations (16) form a six-dimensional system, and in terms of the eigenfre- 



quencies u the characteristic equation takes the form [9] 



uj 2 — 7/7 + a/3 + k 



0, 



(19) 



showing explicitly the existence of zero-modes even upon restriction of the problem to the 
equatorial plane. 



Choosing initial conditions t(0) = f(0) = 0, the periodic solutions (15) satisfy 



-an„ 



-rjn c 



Xfi — a/3 — k 



M 1 



6A/ 
R 



it: 3 1 



3M 
R 



(20) 



whilst n* = = n r s = 0. As the period of the geodesic deviation (15) -which can be 
interpreted as the relativistic generalization of an epicycle- differs from that of the circular 
orbit we started from, the point of closest approach (the periastron) shifts during each orbit 
by a fixed amount. This accounts for the well-known precession of the periastron in general 
relativity [H [9] . 

The zero-modes correspond to secular motions described by linear functions 



IK 



(21) 



with 



0, 



(22) 



Observe, that a non-zero value for v r would have been unacceptable, as it would contradict the 
boundedness of the orbits described. However, non-zero values for v l and v v do not cause such 
problems and are perfectly allowed. In fact, such solutions are required for at least two reasons. 
First, as observed in [9], the inhomogeneous terms in the higher-order geodesic perturbation 
equations generate Poincare resonances which have to be removed by such secular terms. 
This is discussed in detail below. Second, even in the first-order approximation the periodic 
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solutions suffer from the problem that the angle and time between periastra can not be 
matched correctly for eccentric orbits; this mismatch accumulates and grows without bound 



in due course of time, unless corrected by the secular terms (21), (22). 

Whilst a non-zero value of A^ is required by the position and precession of the periastron 
and implies non-zero values for (v t ,v tf> ), the shifts in the cyclic co-ordinates A^ and A^ only 
change the origin of time and azimuth angle and are therefore arbitrary. In view of our initial 
conditions we take A^ = A£ = 0. As A r n then is the only remaining relevant component, 



from here on we will for simplicity write A^ = A n . The secular solutions (21) are further 
restricted by the normalization of the four-velocity. To first order this restriction takes the 
form 

u, L — = => eo^-4^ = 0. (23) 



Taking into account eq. (22) the solutions for the secular velocity terms are 



v l = —. A n , v v = —± A n . 24 

ap — 7/7 ap — 7/7 

Furthermore observe, that application of the normalization condition Q implies 




7/(0) = n r per (0) + n r sec (0) = < + A n = y 1 - — . (25) 

Combining all results, the solutions for the first-order perturbed geodesies describing bound 
motion become 

r a OLKT 

t = — , an r srnwr H crA n , 

3M uj ap — 77/ 

R 

R + an r c cos wr + crA n , ^26) 

V r ■ . '/'- a 

<p = \ 1 = crn c sinwT H crA n . 




uj aj3 — 77/ 

We now relate the parameters A n and a to observable quantities. First, note that the periastra 
and apastra of the orbit occur at proper times r n such that ut„ = nn. We take the even 
values of n to correspond to closest approach (periastron, pa) and the odd values of n to 
maximal distance (apastron, ad). Then we have 

r pa = R + a(A n + n r c ), r aa = R + a (A„ - ri c ) . (27) 

Inverting these equations we get 



In view of eq. ( 25 ) it follows that a takes the value 



hence the dimensionless parameter a/R is related to the eccentricity of the orbit. Starting 
from these expressions we can calculate the energy and angular momentum per unit of mass 
for these perturbed orbits: 

= + (l-^)Jp, i n = i + 5i = r 2 ^, (30) 

to find the first-order change in the constants of motion as compared to the circular orbits: 

1 — = = -uj 2 RA„. (31) 

It follows directly, that the values of (e, £) are unchanged if and only if the secular contributions 
vanish: A n = v = = 0. However, in general we want to allow changes of these values, as 
this is required to describe also non-circular orbits. 

More precisely: up to orientation bound geodesies in Schwarzschild space-time are char- 
acterized by the two constants of motion e and I. Together these two parameters determine 
the angle as well as the time lapse between successive periastra. However, for circular orbits 
these parameters are not independent, as both are determined completely by the value of 
the radial co-ordinate R. Therefore it is not possible to choose a zeroth-order geodesic which 
is circular and has arbitrary independent preassigned values of e and I. In contrast, at first 
order it becomes possible to adjust the parameters such that e and I can be given independent 
values, but only by changing these constants of motion compared to the values they have on 
the circular geodesic. Therefore in general one has to choose a non-zero value for the secular 
contributions from non- vanishing A n . In practice it is easiest to fix the periastron parame- 
ters (radial distance, angle and time lapse between successive periastra) to have preassigned 
values, relating them to e and I afterwards by the procedure described above. 

4 Second-order perturbations 

The first-order geodesic deviations are solutions of the coupled homogeneous linear differential 



equations (14) and combine both periodic and secular terms. In the second-order deviation 
equations these solutions reappear in the inhomogeneous terms. This is evident from the 
second equation ([5]), which can again be cast in a non-covariant but more tractable form in 



terms of = — T^n^n 1 



+ 2u*f AI T + u^d u TJ m» = S» [n] , 

(32) 

dr dr ° K dr 



S"[n] = -2f Al T — — - 4d K f J u x n«— - d a d K f^u\^n 



As expected from eq. ^ the left-hand side of this equation is identical to the first-order 
equation (14). Thus an arbitrary solution of eq. (14) can be added to any solution of (32). 
The right-hand side of these equations is a quadratic expression in the first-order solutions 
and their derivatives. As in general the first-order solution is a combination of periodic and 
secular terms, the inhomogeneous terms S^[n] will contain various products of periodic and 
secular terms. For example, considering again the case of a reference geodesic along which 
the metric is constant, these inhomogeneous terms will be of the form 

5 M [n] = A£ cos 2ut + sin 2wr + B% cos lot + B% sin wr + C\ (33) 
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where oo is an eigenfrequency of the linear operator acting on m M on the left-hand side of the 



first equation (32). When the coefficients B£ s are non-zero we have resonant driving forces 
which lead to singular results for the amplitude of . 

The origin and resolution of this kind of singular behaviour in the perturbative treatment 
of non-linear oscillators was realized long ago by Lindstedt and Poincare [T2] (for a modern 
presentation, see e.g. [2]). Briefly, the dependence of the frequency on the amplitude of an 
anharmonic oscillator is not properly taken into account by the naive perturbative treatment. 
An improved perturbation theory can be developed in which both the amplitude and the 
frequencies of the perturbative solutions are made to depend on the expansion parameters, 
so as to cancel singular behaviour of the final solutions. 

In the present case the procedure is a little more involved, as we have to solve perturbative 
equations of motion (geodesies) subject to a constraint: the normalization condition u 2 = — 1. 
Therefore we proceed as follows. We start from the given reference geodesic x^(r), its tangent 
vector u^t) and the metric on the reference geodesic g^ v = g^ u {x) and its derivatives as 
given. Then we develop the series expansion of neighboring geodesies x^t;^"] as in eq. ([2]); 
both curves x m (t) and x^{t) being geodesies, substitution of this expression into eq. ([!]) gives 



D 2 x fl 



d 2 n^ x - u dn v x „„ - u ,. 



n 



(34) 



+ 



Now define the new time variable A by 



ooX 



00T 



oo + au>± + a 0J2 



(35) 



where oo is the eigenfrequency characterizing the first-order deviation, and wi,2,... are higher- 
order corrections. The expansion then takes the equivalent form 







a 



d?n^ 



,dn v 



+ 2° 



d 2 m^ \- ,,dm u 

-dV +2uT ^la 



(36) 



m 



n 



+ ..., 



where the inhomogeneous source term for is changed to 



[n] 



= ^ dn x dn u 
Xu dX dX 



Xu 



u X n K 



~dX 



^„,X„,U <7„« 



u u n n 



Auii d 2 n^ 4wi - „ ^dn v 
i if ^ 7/ A 

oo dX 2 oo Xu dX ' 



(37) 



It is now possible to choose oo\ and its higher-order generalizations so as to cancel the dan- 
gerous contributions in the inhomogeneous terms S^[n] which produce the divergences. Ob- 
serve, that this is achieved by a rearrangement of the perturbative expansion of the geodesies 
2^[t;o"]. In the next section we illustrate the procedure for the example of bound geodesies 
in Schwarzschild space-time. 



7 



5 Second-order geodesies in Schwarzschild space-time 

In this section we construct the solutions to the second-order deviation equations for bound 
orbits in Schwarzschild space-time. However, we must first briefly return to the first-order 
deviations. The equation for the first-order deviations is the linear homogeneous equation 



dX 2 + 2u \v dX +u « u kX 



0., 



with regular periodic solutions 



n perW = n c cos wA + sinwA. 



(38) 



(39) 



Substitution of this expression in eq. (38) returns eqs. (16), showing that the first-order 
geodesic deviations of a circular orbit are unchanged in the rearranged perturbation theory, 
except for changing r — > A: 



n*(A) 
n r (X) 



a anX 

- — n c sin loX H A r 

to ap — 77/ 



rf c coswA + A n , 



— n r sin a; A + 



(40) 



—a A n, 

ap — 777 



with ujX = cut, showing that the proper frequency has been changed by terms of order a and 
higher. Also note the relation 



V«V = -eon 1 + eon* = 0. 



It follows directly, that n^(A) is an exact solution of the condition (23) 
„Dn v 



9^ 



Dt 



.. ( dX dn u x - ,. „ 

^ u \Trlix +uT >> 



_ Dn» dX 
9 ^J)Xdr= > 



where the pull back of the covariant derivative is 



Dn^ d-nl 1 dx } 
+ 



DX 



dX dX 



Xv 



n 



dn 11 dr 
~dX + dX 



-u T^n . 



(41) 



(42) 



(43) 



Substitution of the solutions (40) in the inhomogeneous terms (37) of the second-order devi- 
ation equations now provides expressions for of the form 

E* = a* sin2a;A + b l sinwA, 

a r cos 2uj X + V cos uj X + c r , (44) 
a v sin 2ojX + b v sin coX. 



TP 



The coefficients (a^,b^,c^) are complicated expressions in terms of R and the first-order 
deviation parameters (n£,A n ), which are given in appendix [B] Then the solutions for 
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take much the same form, with additional secular terms 

m* = m\ sin 2a; A + m\ sin a; A + w t X, 



m r = m r 2 cos 2u)\ + m\ cos a; A + A m , 

mf = m 2 sin 2a; A + mf sin a; A + w^X, 
where the coefficients are solutions of the linear systems 



— 2coa 
-(4a; 2 + k) 
—2a; 77 



(45) 






( m\ \ 






) 


rrf 2 








V m 2 J 




\ <*■ ) 



—ujo. 
(a; 2 + k) 

—LOT] 



-u 2 



m\ 
ml 




and 



(3w l - jw v 



(46) 



(47) 



(48) 



The subtlety in solving these equations is, that the determinant of the matrix of coefficients 
of the nil vanishes as a result of the relation ( 19 ). The system of equations (47) can be solved 
only, if the inhomogeneous terms W have vanishing components in the direction of the zero 
mode of the coefficient matrix; for this to happen, the frequency shift uj\ has to be chosen 
properly. 

To be precise, for u 7^ the coefficient matrix on the left-hand side of eq. (47) has three 
different eigenvalues: (0, — a> 2 , —(2a; 2 + n)). The left zero mode is given (up to an irrelevant 
normalization factor) by 

m M = (P,u, -7)- (49) 



Therefore the condition that eq. (47) is invertible is 

+ ub r - ~fW = 0. 



Now all coefficients are of the form 

If 



n„ 



An 

R 



(50) 



(51) 



with (F^jG 1 *) as given in appendix [Bj determined only by M and R. As a result we finally 
obtain 

oj PG t + coG r - jG^ R ' [ ' 

It should be noted, that the frequency shift u)\ is proportional to the secular radial shift A n , 
but independent of n£, and hence vanishes whenever A n = 0. Provided the constraint (50) 
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on the inhomogeneous terms is satisfied, eq. (47) can be inverted to yield 

rjb 1 



m\ 



m- 



■7 + 



- 7^ 



uj 2 (oj 2 + k) (2u 2 + k)(uj 2 + k 

- -yW 
u{2uo 2 + kY 



at, 



(53) 



r?6* - aW 
uj 2 (u} 2 + k) 



(2uj 2 + k)(uj 2 + k) 



By construction these coefficients satisfy the same constraint as the source terms b^: 

f3m\ + tom\ — jmf = 0, 



(54) 



i.e. they are also orthogonal to the zero-mode (49). In contrast, eqs. (46) for m% can be 
inverted straightforwardly, with the result 



m\ 



m 



rn o 



a 



i2w4 {pJ + W 
-^3 (f3a t + 2ua r 



7a 1 



7a* 



7a Y 



a 

4^2' 



(55) 



4w 2 ' 



Of course, these solutions are determined only up to a solution of the homogeneous 
equation, depending on the initial conditions. 

As concerns the normalization of the four- velocity, in stead of eq. ( 23 ) the solutions up to 
second order must now satisfy the condition 



a u 



+ 



Dn» Dn v 

15715^ 



0. 



(56) 



where = m fI + T^n^n". We have already seen in eq. (42) that the first term vanishes 
identically if we substitute the solution (40) for n M (A). Upon using the first equation ^ the 
term of second order in a in (56) vanishes if 



Dk>* 1 d 2 n 2 
+ 



" Dt ' 2 dr 2 



dr T M + 2 dr 



0, 



(57) 



where n 2 = g^ u n^n ly ' . But to this fixed order in a we are free to replace the proper time by 
the evolution parameter A. After substitution of the explicit expressions for u^, and fc M 
eq. (57) then reduces to 

, dm? 



£0 



dm 



dX 



2lo (eom\ ~~ ^o m 2^) cos ^ A + uj [eom\ — iomf^j cos a; A + Eqw 1 — Iqw 1 * 



oo 2 {n r c ) 2 o x 2oj z A n n r c x 3 u 2 (A 

COS ZOJ A 577- COSUM 



2M 
R 



2M 
R 



4 1 



3M ' 
R 



(58) 
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The two relations that follow for m^ 2 by comparing the terms proportional to coswA and 
cos2wA are identities, implied by eqs. (53) and (55); however, for the constant terms we find 
a constraint 



Eqw* - £ w v 



3 0J 2 (A n y 



4 1 



3ili 
R 



(59) 



Together with the relation (48) this can be used to express and in terms of A m and 
the lower-order parameters: 



w 



3 oj\A n ) 



1 R(c r + K A m ) 



3/2 



3M \ 
R 



3 M w 2 (A r 



AR\ R 



3M 



(60) 



1 R {c r + nA r 



1 



3iU 



3/2 2 V M 



3Ai 
R 



with c r as given in eq. (Ill), and A m a new free parameter to be fixed by the boundary 
conditions on r. More specifically, the net result for the geodesic solutions to second order is 

t = UqT + U{ sin cot + U\ sin 2u)t, 

r = £7£ + C/fcoswr + t/ 2 r cos 2cjt, (61) 

cp = UqT + Uf sinwr + XJ^ sin 2Qt, 
where to this order ui = uo + au)i, and 

I aaKLuA n a 2 uiw t „+ <rcm!! 1 



+ 



yJl-^M. u(aP - 7??) 
1 



+ 



2w ' 



C7? 



cr 2 m 2 , 



U r =R + aA n + ^ a 2 A m , [/[ = ct< + - ct 2 tt^ , 



1 M 



U. 



L 2 r 
2 " g CT m 2' 



ar]KU)A n + cr 2 a;u> y ^ 



crqn c 



1 



^1 _ 3M w(a/3 - 77?) 



2w 



2 V 



(62) 

The second equation (61) implies, that again the periastra and apastra occur at proper times 

(63) 



UJT n = 77 7T. It follows that 



' pa 



t/j - c/[ + u r 2 . 



Then the equations ( 28 ) are modified to 

1 



1 



(r pa - r aa ) = crn r c + - a m{, 



^ {rpa + r aa ) - R = cr A n + i cr 2 (A m + 7773) 



(64) 
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where still relation (25) holds between n r c and A n . Hence these equations determine A n and 
A m in terms of the observables (r pa ,r aa ), and the expansion parameter a: 



a A, 



1 



2 r 

a m l 



1 



' pa 



2M 1 2 /a r L 2M 

— , -cr (A m + +m 2 ) = rpa-R-adl - — . 

(65) 

The parameter a itself is fixed by the time or angle between periastra, which can be expressed 
in terms of the proper-time lapse between periastra: 

UT2 = (oj + (tuji) At = 2tt. (66) 

As in the first-order approximation we can fix the constants of motion e and £ by evaluating 



the expressions (91) after substitution of (61) to second order in a; the result is 



MU%(U{f MuU{U{ 



+ 



(67) 



i m = \u${UZf + \U${U{f + uUr Q U[U*] 2 , 

L J cr z 

where the notation implies that the expressions are to be truncated at order a 2 . 

6 Numerical results 



In the previous sections we have constructed a covariant perturbation theory for orbits of test 
masses in curved space-time based on the method of geodesic deviations, which provides an 
alternative to the common post-newtonian expansions. We have applied it in particular to 
motion in a Schwarzschild background geometry, and obtained explicit expressions for bound 
orbits to second order in the geodesic deviation parameter. Using these expressions for the 
orbits, we can now investigate how well the resulting epicycle orbits compare with the ones 
that are obtained when the geodesic equations are solved purely numerically. As advocated, 
the main advantage of the geodesic deviation method is the fact that the curvature of the 
spacetime is taken into full account, and hence it is expected that the orbits will remain 
accurate even when considered very close to the black hole. In order to test this, two explicit 
examples will be considered of orbits in the region of extremely curved spacetime, i.e. that 
have their periasta close to 6M, the radius of the innermost stable circular orbit (ISCO). 
The comparison to the purely numerically calculated orbits will be done for both the epicycle 
expressions up to first order, and up to second order. 



For the first example the mass of the black hole is set to M = 10 in some unspecified 
units. In the same units, we take the following periastron distance r pa , apastron distance r aa , 
periastron shift 5ip and proper timeshift At between successive periastra: 

r pa = 63.1228, r ap = 92.5279, 5ip = 7.92000, At = 2984.06 (68) 

(It should be noted that the fact that this orbit has a periastron shift greater than the 
Newtonian angular distance between periastra of 2ir, is evidence of the very strong relativistic 
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effect this close to the central mass.) This orbit is uniquely fixed by the following values of 
the energy and angular momentum per unit mass of the test particle: 



0.94827 



35.5000. 



(69) 



This orbit will now be constructed up to first and up to second order in epicycle perturbation 
theory. In the case of the first-order orbits, eq. (26), there are three constants that need 



to be assigned values: a, R, A n . This means that these functions can be subjected to three 
boundary conditions to fix the orbit. The following will be used: the orbit must have a 
periastron shift S(p, have two successive periastra a proper time At apart, and yield a radial 
periastron distance r pa . Prom eq. (26), these three conditions mathematically translate to 
the algebraic conditions 



6<p 
At 




+ 



■y _ 3M uj aft — 777 
R 



aA r 



R + a{n r c + A n ) 



(70) 



Solving the three boundary conditions then yields the following values for the normalisation 
parameter a and the independent epicycle parameters R, A n as well as the associated value 



of n r c : 



a 



-25.9415, 



R = 85.8422, 



A n = 0.449635, 



0.426159, 



(71) 



and hence the following expressions for the orbital functions t(T),r(T),ip(T): 

t(r) = 1.2851 r + 17.9319 simW) 
r(r) = 74.1780- 11.0552 cos(wr) 
cp(r) = 0.00611 r + 0.46944 simW) 



in which the epicycle frequency is given by 



UJ 



0.00270. 



(72) 



(73) 



Using the expressions for the e n and £ n from eqs. (30) and (31 ), the energy per unit mass and 



the angular momentum per unit mass of this epicycle orbit are given by 



0.94644, 



35.1840, 



(74) 



which are at most a few tenths of a per cent different from the values in eq. ( 69 ) 



Moving on to the second order, the orbital functions eq. (62) can be subjected to the 
same three boundary conditions as before, but this time an extra constant A m appears that 
needs to be assigned a value. This means that a fourth boundary condition can be imposed, 
which will be chosen to fit the radial apastron distance r aa . The four boundary conditions are 
formulated as algebraic conditions by observing that the periastra and apastra correspond to 
the extreme values of r(r) in eq. (61) and that these occur at times ujT n = mr, n an integer 



number. From <f>{r) then follows the periastron shift, while r(r) yields the radial distances 



13 



of the periastra and apastra. The four boundary conditions then translate to the following 
algebraic conditions: 



UJ + aUJ\ \0J + (TU3\ I 

r pa = U r v+U[ + U r 2 , r aa = U r -U{ + U T 2 . (75) 

Solving these conditions then yields values for the normalisation parameter a and the in- 
dependent epicycle parameters R, A n , A m as well as the associated value of rf c . However, 
when comparing the epicycle approximation in higher order to a specific known orbit, there 
is generically more than one set of values <7, R, A n , A m solving the boundary conditions. The 
question then presents itself which of these orbits is the most accurate one to describe the 
(unique) solution to the geodesic equation parametrized by the values e, I. The answer is 
provided by the previously made observation that the expansion parameter of the epicycle 
approximation is the dimensionless number an r c /R, and hence that the solution set that pro- 
duces the smallest value of this parameter yields the most accurate orbit. Using this criterium, 
the following solution set a, R, A n , A m and associated value n r c is chosen: 

a = -3.11032, R = 79.8763, A n = -2.88620, A m = -2.19232, n r c = 3.75194, (76) 

which give rise to the following expressions for the orbital functions 

t(r) = 1.28429 r + 22.0076 sin((w + fjwi)r) + 1.47690 sin(2(w + CTWi)r), 
r(r) = 78.2490 - 14.7025 cos((w + <twi)t) - 0.42359 cos(2(w + aui)r), 
<p(r) = 0.00611 r + 0.57279 sin((w + fjwi)r) + 0.04272 sin(2(w + CTWi)r), 

in which the frequencies are given by 

u = 0.00279, erwi = -0.00009, u + ouji = 0.00270. (77) 



From eq. (67), the energy per unit mass and angular momentum per unit mass for this 



approximation to the orbit then follow as 

e m = 0.94807, l m = 35.5802, (78) 



which are at most a few hundredths of a per cent different from the values in eq.(69) 



Having found the epicycle expressions to first and second order, they can be compared to 
orbital functions (t(r), r(r), <^(r)) as calculated by solving the geodesic equations by purely 



numerical means. The latter are completely specified by the values e and I of eq. (69) by 
methods explained in |16j . 

In Figure [TJ the radial function r(r) up to first and second order is given divided by its 
purely numerical counterpart, as a function of proper time. As can be seen, the relative 
difference between our first-order approximation and the numerical one is at most about 
8%; introducing the second-order epicycle improves the relative difference to less than 0.5%. 
Figure [2] shows the absolute difference between the angular co-ordinate tp in the epicycle 
approximation (both to first and second order) and in the numerical one. As can be seen, the 
first-order epicycle deviates from the purely numerical one by at most 0.4 radians during any 
period, whereas the second order deviates from the numerical one by at most 0.25 radians. 
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2000 4000 6000 8000 

Figure 1: The radial function r(r) in epicycle approximation up to first (solid line) and second 
(dashed line) order divided by the numerical one, as a function of proper time r. This is the 
result for the first example mentioned in the text. 

0.6 
0.4 
0.2 

o.o 

-0.2 
-0.4 
-0.6 

Figure 2: The difference of the angular function ip(r) in epicycle approximation up to first 
(solid line) and second (dashed line) order minus the numerical one, as a function of the 
numerical ip in units of 7r.This is the result for the first example mentioned in the text. 

The figures show that the radial coordinate r(r) is very well approximated by the second- 
order epicycle functions. On the other hand, in the second-order epicycle approximation the 
coordinates <p(r) and t(r) still show some small deviations; further improvement could be 
made by including third- and higher-order corrections. 

As a second and less extreme example, the mass of the black hole will again be set 
to M = 10, and an orbit will be considered which has, in the same units, the following 
periastron distance r pa , apastron distance r aa , periastron shift 5ip and proper timeshift At 
between successive periastra: 

r ap = 100.000, r aa = 150.000, 6<p = 2.61000, At = 3373.56. (79) 

This orbit is bound between 10M and 15M, and corresponds to one that has energy per unit 
mass and angular momentum per unit mass given by: 

e = 0.96362, I = 40.0892. (80) 
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Subjecting the first order epicycles of eq. (26) to boundary conditions similar to the ones 
discussed in the previous example, the values for the normalisation parameter a and the 
epicycle parameters R, A n as well as the associated value n r c follow as: 

R = 125.515, a = -27.8282, A n = 0.06148, n r c = 0.85540, (81) 

which give the following first order orbital functions: 

t(r) = 1.14879 r + 22.1255 simW), 
r(r) = 123.804- 23.8041 cos(o;r), 

(p(r) = 0.00264t + 0.52501 sin(wr), (82) 
in which the epicycle frequency is given by 

U = 0.00186. (83) 



Using the expressions for the e n and l n from eq. (31), the energy per unit mass and angular 
momentum per unit mass of this epicycle orbit are given by 

e n = 0.96325, £ n = 40.4225, (84) 



which deviate from the numerical counterparts of eq. (80) by less than a per cent. 

Going to second order eq. (62), an analysis similar to the one discussed in the previ- 
ous example yields the following values for the normalisation parameter a and the epicycle 
parameters R, A n , A m as well as the associated value for n T c : 

a = -20.6911, R = 120.507, A n = -0.23446, A m = 0.00745, n r c = 1.14772, (85) 

which give the following second order orbital functions: 

t(r) = 1.14876 r + 23.2580 sin((w + fjwi)r) + 2.31954 sin(2(w + CTWi)r), 
r(r) = 126.953 - 25.0000 cos((u + <rwi)r) - 1.95325 cos(2(w + aui)r), 
<p(r) = 0.00264 r + 0.55205 sin((w + fjwi)r) + 0.06398 sin(2(w + CTWi)r), 

(86) 

in which the frequencies are given by 

u = 0.00195, <rui = -0.00009, u + cr^i = 0.00186, (87) 
The energy per unit mass and the angular momentum per unit mass follow from eq. (67) as 

e m = 0.96362, £ m = 40.0822, 



which deviate less than a hundredth of a per cent from their purely numerical counterparts 
of eq. ((80|). 

Having found the epicycle expressions to first and second order, they can again be com- 
pared to orbital functions as calculated by solving the geodesic equations by purely numerical 
means. The latter are completely specified by the values e and I of eq. (80) by methods 
explained in [16], the result being numerical functions r(r) and <p(t). 
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In Figure [3j the radial function r(r) is shown up to first and second order as a function of 
proper time, compared to the purely numerical approximation in terms of the ratio. As can 
be seen, the relative difference between our first-order approximation and the numerical one 
is at most about 4%; introducing the second-order epicycle improves the relative difference 
to less than 0.5%. Figure [4] shows the absolute difference between the orbital function ip 
(both to first and second order) and the numerical one, as a function of the numerical (p. 
As can be seen, the first-order epicycle deviates from the purely numerical one by at most 
0.1 radians, whilst the second-order one deviates from the numerical one by at most 0.02 
radians. These results imply that the orbital functions as given by the geodesic deviation 
method very accurately describe the geodesic orbit uniquely specified by the values of e, t as 
given in eq. ( 80 ) . 




Figure 3: The radial function r{j) in epicycle approximation up to first (solid line) and second 
(dashed line) order divided by the numerical one, as a function of proper time r. This is the 
result for the second example mentioned in the text. 
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0.00 r 



-0.05 
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Figure 4: The difference of the angular function <p(t) in epicycle approximation up to first 
(solid line) and second (dashed line) order minus the numerical one, as a function of the 
numerical ip in units of 7r.This is the result for the second example mentioned in the text. 
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7 Conclusions and outlook 



The two examples presented show that the geodesic deviation method can be used to ac- 
curately describe eccentric bounded geodesies of a Schwarzschild black hole. The choice of 
boundary conditions used in this article focussed mostly on improving the accuracy of the 
radial function and, as such, the method effortlessly produces very accurate results for r(r), 
even when describing orbits grazing the ISCO. The accuracy of t(r) and <p(t) can be further 
improved by taking the method to higher order, or for specific purposes by choosing boundary 
conditions that focus on constraining these functions rather than the radial values. Although 
each extra order will introduce an extra degree of freedom and hence allows for one extra 
boundary condition to be used as a constraint on t(r) and <p(t), it is expected that only one 
extra order correction suffices to improve both two functions in accuracy. This is because of 
the normalization of four- velocity and the fact that the radial orbital function is already very 



accurate at second order: the relation eq. (90) 



dr J V r ) V dr J \dr 

will ensure an increase in accuracy of t(r) or (p(r) when the extra degree of freedom is used 
to improve the accuracy of the other. 

At least two further applications of the geodesic method are immediately foreseeable. 
Firstly, it opens up the possibility of removing one numerical step in calculating the grav- 
itational radiation in a Schwarzschild background. A formalism to calculate gravitational 
radiation due to the motion of a test mass in a Schwarzschild background has been topic of 
research for many decades and was presented in a final form in [T7] , [18] , but where that work 
had to rely on numerical descriptions of the geodesic orbits, the geodesic deviation method 
allows to replace that numerical step by an analytical one. Secondly, our presentation of the 
geodesic deviation method only relied on the assumption that the equations of motion of the 
deviation vectors were linear differential equations with constant coefficients, and the results 
are therefore quite general and can be used to describe other physical situations as well. For 
example, the method can be applied to the description of the motion of an electrically charged 
test mass in a Schwarzschild spacetime while experiencing a Lorentz force due to an axially 
symmetric magnetic field. Both of these applications will be the topic of future publications. 
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A Geodesies in Schwarzschild space-time 



In this appendix we collect some well-known result about Schwarzschild geometry which are 
relevant for the results derived and discussed in the main text. We describe the static and 
spherically symmetric exterior geometry of Schwarzschild space-time in terms of the standard 
Schwarzschild-Droste co-ordinates (in units in which c = G = 1): 

dr 2 =[l- 2 ^ dt 2 - - r 2 (d9 2 + sin 2 9 dip 2 ) . (90) 

Considering geodesies in the equatorial plane = 7r/2, there are two constants of motion 
representing the energy and angular momentum per unit of mass: 

e=(l-™)* i = r 2 ^. (91) 
V r ) dr dr 



Using eq. ( 90 ) the radial motion is then described by the equations 

dr\ 2 _ 2 / 2M \(i £2 \ d2r - M £2 (i 6M \ 
dr J \ r J \ r 2 ) ' dr 2 r 2 r 3 \ r ) 



(92) 



It is not possible to solve the general equation for r(r), except in special cases like circular or- 
bits. However, one can exchange the proper time dependence for azimuth angular dependence 
by using the second equation (91), and determine r((p) from 

dl Y-e 2 (l 2M \(l+ l2 \ (93) 
dip r J \ r J \ r 2 I 



I 2 



This equation has solutions in terms of elliptic integrals. On the other hand, the geodesic 
deviation method allows to find approximate analytic solution for r(r) itself starting from a 
special solvable geodesic, such as the circular orbit. 

Even though it does not produce explicit solutions for r(r), the shape equation (93) 
contains useful information. To extract this information, we introduce a new variable y |15j 
such that 

r=— ~ , (94) 

1 + e cos y 

where the parameters (e, a) are implicitly given by 



(95) 



at 2 

a 



1 w +f ( 3+e2 )= - 



In the newtonian approximation e represents the eccentricity of the elliptic orbit, and (1 — e )a 
is the length of the semi-major axis. Eq. (93) implies the full solution to satisfy 

dy\ 2 2M 

— =1 (3 + ecosy). (96) 

dip J a 
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Clearly, for large o>M the second term is negligeable and we can take y = <p — ipo, which 
returns the exact newtonian orbit. Close to the horizon this is no longer applicable, as is shown 
for example by the existence of an innermost stable circular orbit (ISCO) with r = 6M. 



The periastra of the geodesies ( 94 ) are reached for y = 2im. Eq. ( 96 ) shows, that the period 
in y is not identical to the period in <p, explaining the periastron shift in Schwarzschild-Droste 
co-ordinates. We can now derive integral expressions for this periastron shift in azimuth angle 
ip and in observer time t, as follows. The total change Acp between two periastra is 

Aip = 2ir + 5ip = dy 1 (97) 

Jo ^l-2M(3 + ecosy) 

where 5ip is the advance of the periastron compared to the previous one. Now using the 
conservation laws we can also write an expression for the time lapse between two periastra 

ea 2 r2n 



At=— dy . (98) 

I h ( 1 + ecoBy )2^1_^(l + ecosy )j^l_2M( 3 + ecOB y) 

It follows plainly that the angular shift Aip and the time lapse At between two successive pe- 
riastra are independent orbital characteristics, determined by the two independent parameters 
(e,£) or equivalently (a, e). 

The counting of parameters determining geodesies is simple: once geodesies are restricted 
to the equatorial plane = it/ '2, there remain three ordinary second-order differential equa- 
tions for the co-ordinates (t,r,tp). The solutions dependend a priori on six constants of 
integration: three for the co-ordinates and three for the velocities. However, the velocities are 



restricted by the normalization of the four- velocity u 2 = —1, cf. the first eq. (92). Therefore 
in practice we fix the origin of time and azimuth angle, the initial value of radial co-ordinate - 
which we usually take to be the radial position of the periastron- and two velocity parameters 
represented by the constants of motion (e, £). 



B Driving terms for the second-order deviations 

In this section we present the explicit expressions for the coefficients {a^,b^,c^) in the ex- 
pansion of the driving terms S^ 4 for the second order deviations, eq. (44), in the case of 
Schwarzschild geometry. The general expression for these terms is 



^ = 5^ + — T", 



with as given in (32) and (33), after substitution r — > A, and defined by 

. d 2 n^ 



d\ 2 



4r 



dn v 
~d\' 



(99) 



(100) 



In the expansion of S 1 ^ as defined in eq. ( 33 ) : 



S>*[n] = A£ cos 2wr + sin 2wr + B£ cosur + B% sin lot + C, 
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the coefficients A c = A r s = A? = vanish, whilst 



At _ a t _ 1 fn n2 

~~ " ~~ 1 2M 



ii r o M 18A/ 2 

i? 4 (i- (i - m\ ' 



R J 



R 



Next, B l c = B T s = Bf = 0, and 



13M 6M 2 * 



#2 = ^77w T77V 4--— + —=- A n <, 



(l-¥) 



" ^ 4 (i-^) (i-^) 2 V * ^ 2 1 n c 



^ i? i _ 3M ^ # ; A «n c . 

In addition there is a constant term in the r-component: 

, M i , M , u ir 34A/ i 75M 2 54M 3 



The Lindstedt-Poincare correction terms (100) are of the form 

T» = D^smuX + D^ cosuX + E^, 

with L>* = D r s = Df = and 

D l s = — 2wanc, 



whilst E l = = and 



£ r = -2k A n . 



Combining the results to compute 



s > 

OJ U! CO 
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we find the coefficients to be of the form (51) with 

F* = -uja 



13M 
R 



6M 2 



1 



2AI 
R 



(l-^f) 



2M 7 



36M , 48A/ 2 



(l-^f) (l-^) 



2 ' 



-w?7 ■ 



12M 



and 



G 1 
G r 



3M 



2wa, 
4M 



3M 



It then follows by eq. (52), that 



3 1 

2(T- 



i? 3 1 
2urj. 



10M , 18M 2 
-R H 2 



3M \ 
R J 



An 

(!_6M) R 



Note that, whilst on the ISCO R = 6M the fundamental frequency vanishes: to 
Lindstedt-Poincare frequency shift for finite A n is singular there: loi — > oo. 
Finally, the constant c r in the driving terms S r is given by 



C r + — E r . 

UJ 



Using the value (109) for the ratio uj\/uj, this leads to the result 



3M 1 + f . 2 

#4 j_ _ 3M ^ t" 2R2 
H 



3M / A 



R 



26M i 165Af 2 _ 
R T R 2 



:;'!6M 3 _|_ 324M 4 



6A-A 
R J 
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